arXiv: 1506.04310v3 [cond-matsoft] 13Jul2015 


The order-disorder transition in model lipid 
bilayers is a first-order hexatic to liquid phase 
transition 

Shachi Katira ^ Kranthi K. Mandadapu ^ Suriyanarayanan Vaikuntanathan ^ Berend Smit ^ ^ and David Chandler II 

* Department of Bioengineering, University of California, Berkeley, CA 94720, USA,'!'Chemical Sciences Division, Lawrence Berkeley National Laboratory, Berkeley, CA 94720, 
USA,^ Department of Chemistry, University of Chicago, Chicago, IL, 60637, USA,§ Department of Chemical and Biomolecular Engineering, University of California, Berkeley, CA 
94720, USA,^Institut des Sciences et Ingenierie Chimiques, Valais Ecole Polytechnique Federale de Lausanne Rue de I’lndustrie 17, Sion CH-1950 Switzerland, and H Department 
of Chemistry, University of California, Berkeley, Berkeley, CA 94720 USA 


We characterize the order-disorder transition in a model lipid bi¬ 
layer using molecular dynamics simulations. We find that the or¬ 
dered phase appears to be hexatic. In particular, in-plane structures 
possess a finite concentration of 5-7 disclination pairs that diffuse 
throughout the plane of the bilayer, and further, in-plane structures 
exhibit quasi-long-range orientational order and short-range transla¬ 
tional order. In contrast, the disordered phase is liquid. The tran¬ 
sition between the two phases is first order. Specifically, it exhibits 
hysteresis, and coexistence exhibits an interface with capillary scal¬ 
ing. The location of the interface and its spatial fluctuations are 
analyzed with a spatial field constructed from a rotational-invariant 
for local 6-fold orientational order. As a result of finite interfacial 
tension, forces of assembly will necessarily exist between membrane- 
bound solutes that pre-melt the ordered phase. 

lipid bilayers, phase transition, hexatic, line tension 

Abbreviations: DPPC (dipalmitoyl phosphatidylcholine) 


In this and the following paper [1], we detail the exis¬ 
tence of a first-order order-disorder transition in model lipid 
bilayers, and we show how pre-melting of the ordered phase is 
responsible for a powerful membrane-mediated force between 
transmembrane proteins. Pure lipid bilayers primarily exist 
in two phases — an ordered phase often referred to as the ‘gel’ 
phase (L /3 or L^/ phase), and a disordered liquid-like phase (Lq, 
phase) [2]. Here, we use molecular dynamics of the MARTINI 
model [3,4] to investigate the nature of the ordered phase and 
the transition to the disordered phase. 

Prior work with the MARTINI model [3, 5] on small sec¬ 
tions of membrane showed hysteresis when transitioning be¬ 
tween ordered and disordered structures. This form of col¬ 
lective behavior is suggestive of singular behavior, but the 
methods of analysis employed in the earlier works were in¬ 
sufficient to demonstrate the scaling inherent to a first-order 
transition and whether the ordered phase is a crystal or some 
other structure. 

In following up on Refs. [3] and [5], we establish that a sta¬ 
ble interface exists between the ordered and disordered phases 
at coexistence, and that the fluctuations of the interface are 
capillary-like [6,7]. These findings establish that the transition 
between the ordered and disordered membranes is first order, 
which is consistent with experiments [8]. Further, we show 
that the spatial structure of the ordered phase exhibits quasi- 
long-range bond-orientational order and short-range transla¬ 
tional order. These findings indicate that the ordered mem¬ 
brane is a hexatic phase, like those considered in melting of 
two-dimensional systems [9-11], but experiments have not yet 
determined whether ordered phases of membranes have this 
structure. 

A molecular model that reproduces the phase behavior of 
a lipid bilayer is needed to understand lipid behavior around 
solutes such as proteins and cholesterol, and therefore to un- 
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Fig. 1. Order-disorder phase diagram in the tension-temperature, A-T, plane. 
The lateral pressure across the membrane is —A. Points are estimated from 10 in¬ 
dependent heating runs like those illustrated in Fig. 2 for a periodic system with 128 
lipids. Insets are cross sections showing configurations of a bilayer with 3200 lipids 
in the ordered and disordered phases. The heads are colored gray while the tails are 
colored pink. Water particles are omitted for clarity. The hydrophobic thicknesses, 
T)o and T>(\, are the average vertical distances from the first tail particle of the upper 
monolayer to that of the lower monolayer. A macroscopic membrane buckles for all 

A < 0. 

derstand lipid-mediated interactions [1,12]. While atomistic- 


Significance 

Lipid bilayers exist in ordered and disordered states — so-called 
“gel” and “liquid” phases. First-order phase transitions between 
the two have long been known, but whether the ordered phase 
is crystal-like or hexatic has not been known. Here, using large 
scale molecular simulation, we demonstrate a first-order transi¬ 
tion and predict that the ordered phase is hexatic, which can 
provide a basis for understanding mobility and organization of 
proteins in the ordered phase. 
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Fig. 2. Structural measures of different phases as a function of temperature, T. (A) Variation in area per lipid with temperature during heating and cooling shows finite 
jumps and hysteresis. (B) Average local orientational order, {(pi), also shows finite Jumps as a function of temperature while heating and cooling. Magnitudes of heating and 
cooling rates are 3 K/iis. (C) Snapshots from the ordered phase at 279 K (top) and the disordered phase at 315 K (bottom). Both display the tail-ends of each lipid in one 
monolayer of the bilayer. The ordered phase exhibits a hexagonal packing, while the disordered phase exhibits a random arrangement of the tail-end particles. Empty regions 
in the disordered configuration are mostly filled by tail-ends from the other monolayer or by non tail-end particles, none of which are not shown in the snapshot. 


level models [13,14] provide fine-scale structural details, they 
are far too computationally expensive to study the phase 
behavior of a bilayer system with any reasonable degree of 
rigor. To enable access to sufficiently large length scales and 
time scales (admittedly at the cost of reduced detail), coarse¬ 
grained models have been constructed for lipid bilayers and 
studied using molecular simulations [3,4,15-17]. Some of 
these models exhibit a possibly solid-like ordered phase, and 
a liquid-like disordered phase. 

We have chosen the MARTINI model because it success¬ 
fully reproduces equilibrium bilayer properties such as the 
area per lipid and bending modulus [4]. Additionally it has 
been applied to studies of the plasma membrane [18], lipid 
rafts [19], organization of transmembrane proteins [20], and 
membrane tethers [21]. Although we choose dipalmitoyl phos¬ 
phatidylcholine (DPPC) as a model lipid, we expect our re¬ 
sults to hold for many other lipid bilayers or monolayers that 
exhibit ordered and disordered phases. 

Order-disorder transition in a model lipid bilayer. Fig. 1 con¬ 
trasts configurations and shows our estimated phase boundary 
between ordered and disordered phases in the DPPC MAR¬ 
TINI bilayer system. The ordered phase has regular tail pack¬ 
ing compared to the disorganized tail arrangement of the dis¬ 
ordered phase. A consequence of the regular tail packing is 
that hydrophobic thickness of the ordered phase, Vo, is larger 
than that of the disordered phase, Pa- Correspondingly, the 
area per lipid in the ordered phase is smaller than that in the 
disordered phase. 

The tension-temperature, A-T, phase diagram shown in 
Fig. 1 was estimated from heating runs of a small system — 
128 lipids solvated by 2000 water particles. Fig. 2A shows 
the change in area per lipid with temperature while heating 
and cooling a bilayer. There are finite jumps in area per lipid 
as the system transitions between the two phases, suggesting 
a first-order phase transition. Hysteresis occurs because or¬ 
dering from the metastable disordered phase is much slower 
than disordering from the metastable ordered phase. Due to 
the difference in time scales, when contrasting melting and 
freezing from heating and cooling runs, the melting points 
from heating runs provide the more accurate estimates of the 
actual phase boundaries. The phase boundary graphed in 
Fig. 1 shows error estimates based upon the statistics of sev- 
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eral heating runs. Systematic errors due to small system size 
and heating rate have not been estimated. 

Viewing the end particles of all the lipid chains in one of 
the two monolayers provides a convenient visual representa¬ 
tion that distinguishes the two phases. (Which of the two 
monolayers chosen for this rendering is irrelevant.) This rep¬ 
resentation is shown in Fig. 2C for the ordered and disordered 
phases. These tail-end particles appear hexagonally-packed in 
the ordered phase and randomly arranged in the disordered 
phase. 

To quantify this impression, we consider Halperin and Nel¬ 
son’s local rotational-invariant [9,10], 

= exp(6i6>y)| , [1] 

jenn{l) 

where Oij is the angle between an arbitrary axis and a vector 
connecting tail-end particle I to tail-end particle j. Here, the 
summation over j G nn(/) is over the six nearest neighbors of 
particle 1. (See inset in Fig. 2B.) The ensemble average, (0^), 
is unity for a perfect hexagonal packing, and it is zero to the 
extent that hexagonal packing is entirely absent. The varia¬ 
tion in ((/);) with temperature for a DPPC bilayer is shown in 
Fig. 2B. The finite jumps in ((j)i) as a function of temperature, 
and the hysteresis, are again suggestive of a first-order phase 
transition. 

Stable interface between the ordered and disordered phases 
exhibits capillary scaling. To establish whether the first-order¬ 
like behavior described above persists to larger scales and thus 
actually manifests a phase transition, we consider larger sys¬ 
tems and the behavior of the interface that separates the or¬ 
dered and disordered phases. Fig. 3A shows this coexistence 
for a system size of V = 3900 lipids with an interface between 
the two phases. The interface, which spans the membrane, 
is equilibrated in the constant area ensemble. This ensemble 
can maintain an area per lipid that is intermediate between 
two phases and can therefore stabilize an interface if, in fact, 
two distinct phases do exist. At such conditions, a line ten¬ 
sion can then be calculated from the power spectrum of the 
interfacial fluctuations. 

To analyze interfacial fluctuations, we first identify the lo¬ 
cation of the interface at each instant. This location is found 
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with a two-dimensional version of the three-dimensional con¬ 
struction described in Refs. [22,23]. Specifically, the interface 
is here defined as a line in the plane of the bilayer where 
the value of a coarse-grained field characterizing local orien¬ 
tational order is intermediate between those of the ordered 
and disordered phases. The field we use is 

= [ 2 ] 

I 

where is the position of the lih tail-end particle projected 
onto a plane parallel to that of the bilayer, and /(r — ri,^) 
is a coarse-grained delta-like function in the two-dimensional 
space, 

/(r - r;;e = (l/2<^) exp (-|r - ry/2^^) . [3] 

Again, the tail-end particles are from the chains of only one of 
the two monolayers forming the bilayer. The field variable, r, 
is a two-dimensional vector specifying a position in the plane 
of the bilayer. The coarse-graining width, is chosen to be 
the average separation between tail-end particles I and j when 
((0i - - {4'i)Y) in the ordered phase is 

1/10. This choice yields a value of ^ = 1.5 nm. 

For numerics, a square lattice tiles the average plane of 
the bilayer, and the coarse-grained field 0(r) is evaluated at 
each lattice node. For convenience, the coarse-graining func¬ 
tion is truncated and shifted to zero at 3^. The instantaneous 
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Fig. 3. Capillary fluctuations of the order-disorder interface: (A) Snapshot of a 
system showing coexistence between the ordered and disordered phases. The gray 
contour line indicates the location of the interface separating the ordered and disor¬ 
dered regions. The snapshot is a top view of the bilayer showing the tail-end particles 
of each lipid in one monolayer. (B) Fourier spectrum of the fluctuations of the in¬ 
stantaneous order-disorder interface. The line is the expected small-A: scaling from 
capillarity theory. 


order-disorder interface is identified by interpolating between 
these adjacent lattice nodes to find the set of py)ints s satisfy¬ 
ing 0(s,t) = (0d + 0o)/2. Here 0d and 0o are (0(r)) evaluated 
in the disordered and ordered phases, respectively. 

Fig. 3A shows a snapshot of the instantaneous interface 
identified in this way. This free interface is stable for the 
length of simulations, l.S fis. For the thermodynamic state 
considered in that case, zero lateral pressure and 294 K, we 
have 0d = 0.4 ± 0.02 nm“^ and 0o = 2.15 ± 0.2 nm“^. 

Fig. 3B shows the Fourier spectrum of the height fluctua¬ 
tions of this interface, {\Shk\‘^). Two different system sizes are 
studied, with the larger system having approximately dou¬ 
ble the interface length of the smaller system. The Fourier 
component Shk is related to the height fluctuation 6h{x) as 
Sh{x) = 5]]^ Shk ex.p{ikx) where x is a point on the hori¬ 
zontal axis in Fig. 3A. Here, 0 ^ x ^ L, and L is the 
box size. With periodic boundary conditions, k = 27rm/L, 
m = 0,±1,±2,---. According to capillarity theory [7, 24], 
~ k^TjL^k^ for small k, with being the Boltz¬ 
mann’s constant. Given the proportionality with 1/A;^ at 
small k (i.e., wavelengths larger than 10nm), comparison 
of the proportionality constants from simulation and capil¬ 
larity theory determines the interfacial line tension, yielding 
7 = 11.5 ± 0.46 pN. (This value is significantly larger than the 
prior estimate of line tension for this model, 3±2 pN [3]. That 
prior estimate was obtained from simulations of coarsening of 
the ordered phase.) 

The stability of the interface and the quantitative con¬ 
sistency with capillary scaling provide our evidence for the 
order-disorder transition being a first-order transition in the 
model we have simulated. 


The ordered phase is a hexatic phase. The ordered phase 
shows hexagonal packing with a large number of unbound dis¬ 
locations, each composed of 5-7 disclination pairs [9]. These 
pairs are shown in Fig. 4A, which renders a 25 x 25 nm^ 
bilayer in the ordered phase simulated at 294 K and zero lat¬ 
eral pressure. By showing three different conhgurations of 
the bilayer separated by 600 ps, we illustrate how the discli¬ 
nation pairs diffuse freely throughout the system. Given that 
a lipid bilayer is a quasi-two-dimensional system, it is pos¬ 
sible that the ordered phase with its unbound dislocations 
could be in the so-called hexatic phase, an intermediate phase 
between the solid and liquid phases in two dimensional sys¬ 
tems [9,10,25-29]. The hexatic phase was first predicted for 
a purely two-dimensional system [10], and it has been pre¬ 
dicted to occur in three-dimensional systems [30]. It has been 
observed in three-dimensional systems, e.g.. Ref. [31], and in 
quasi-two-dimensional systems, e.g.. Refs. [11,32-34]. This 
phase is characterized by short-range translational order and 
quasi-long-range orientational order. 

To test whether the ordered phase of the chosen model 
lipid bilayer system is hexatic, we calculate appropriate trans¬ 
lational and orientational correlation functions. Because any 
finite sample is anisotropic in the presence of quasi-long- 
range order, the correlation functions must be applicable 
when different system orientations are not equivalent. At 
the same time, they cannot be based upon a physical refer¬ 
ence lattice because the system, we shall see, lacks long-range 
translational order and is thus not crystalline. Bernard and 
Krauth [35] showed how to proceed under such circumstances, 
and we follow their approach. 

To do so for our membrane system, we again project po¬ 
sitions of tail-end particles from one of the monolayers onto 
the average plane of the bilayer. Translational order is then 
examined by considering the pair correlation function along a 
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specific sample orientation. The pair correlation function is 

S(r) = +rj)), [4] 

where ( 5 (r) is the delta function for r in the plane of the bi¬ 
layer, and p is the number of particles per unit area. Its eval¬ 
uation along a specific direction, g{x)^ is ^(r) with r = xx, 
where x is the unit vector of an axis of the Cartesian coor¬ 
dinate system for which T = (1/^) '^i is maximal. Here, 
= (l/ 6 )E, 6 nn (0 exp( 6 z^j 7 ), with the summation taken 
over the six nearest neighbors j of particle 1. This pair corre¬ 
lation function, g{x), is denoted by ^(x, 0) in Ref. [35]. 

For a solid with long-range translational order, g(x) will 
exhibit a power law decay, specifically g{x) — 1 ~ x“^, with 
^ ^ 1/3 [36], but our results in Fig. 4B show a more rapid 
decay, seemingly exponential. This rapid decay is indicative 
of short-range translational order characteristic of a hexatic 
phase [9]. Our calculations for the ordered phase shown in 
Fig. 4 are for a system of 51,200 lipids (solvated by 800,000 
water particles) at temperatures T = 294 K. The initial con¬ 
figuration was created by periodically replicating a smaller 
system. By construction, therefore, the initial configuration is 
effectively crystalline, with a large unit cell (or tile). The equi¬ 
librated system with short-range translational order emerges 
only after several microseconds of simulation. See Methods 
and Supporting Information (SI). 

For orientational correlations, we consider 

= exp(6i6';j)<5(r - Fi). [5] 

I jenn{l) 

and compute 

G6(r) = (■!/>6(r)'iA6(0)) . [6] 

We then determine Geix)^ which is ^0(1*) evaluated at r = x x, 
where again x is the unit vector along an axis of the Cartesian 
coordinate system that maximizes T. For the hexatic phase, 
this correlation function will be quasi-long-range (power-law 
decay), progressing to zero more slowly than x“^/^ [25,36]. 
Our results, shown in Fig. 4C and D, are consistent with that 
behavior. While the system size is large, it is not large enough 
to determine the exponent for G 6 (a^)’s decay. Nevertheless, it 
is clear from the figure that orientational correlations in the or¬ 
dered phase are not long ranged, but rather quasi-long ranged 
over the distance scales accessible to our simulations. Further, 
the orientational correlations decay much slower than trans¬ 
lational correlations. This striking difference between trans¬ 
lational and orientational correlations does not exist in the 
disordered liquid phase. In particular, our results in Fig. 4C 
for T = 312 K, which is above the order-disorder transition 
temperature, show short-range orientational correlations. 

The dislocations highlighted in Fig. 4A and the correla¬ 
tion functions in Fig. 4B, C and D are calculated for the plane 
consisting of tail-end particles of one monolayer. We have also 
studied the structure and the correlation functions for planes 
corresponding to each particle along the length of the tail. As 
could be anticipated [37, 38], we observe that the number of 
dislocations decreases as we move upwards in the lipid chain. 
However, despite the disorder gradient, the translational and 
orientational correlations at every level along the lipid chains 
are consistent with the hexatic phase (see SI). The relative 
stability of the hexatic phase over a solid phase in this sys¬ 
tem can be interpreted in terms of the presence of thermally 
induced fluctuations in the midplane of the membrane [39]. 
These fluctuations can cause unbinding of dislocations, which 
destroys translational order. 

Taken together — the finite concentration of diffusing 5-7 
disclination pairs, the quasi-long-range orientational order. 


and the short-range translational order — indicate that the 
ordered phase of the model lipid bilayer we study is a hexatic 
phase. 


Model dependence and relation to experiments. X-ray 
diffraction experiments have shown that an ordered hexatic 
phase exists in between the solid and disordered phases of 
lipid monolayers [32]. In the case of lipid bilayers, the X- 
ray scattering data found a finite in-plane correlation length 
in the ordered phase, about 20 nm [40], but the experiments 
were unable to distinguish between a hexatic phase and a solid 
phase with finite sized domains. Perhaps our prediction will 
motivate further experimental work to resolve this issue. 

Bear in mind that while the MARTINI model is able to 
reproduce the order of the transition in DPPC bilayers, its 
ordered phase does not exhibit the experimentally observed 
tilt characteristic of this phase, and tilt can be induced in 
the MARTINI model by decreasing the size of tail particles 
relative to head particles [3]. Tilt can also be induced in a 
soft-sphere model by increasing repulsion between heads [41]. 
Such features and the coupling between tilt and the order- 
disorder transition can be explored in future work. 

Also bear in mind that our analysis of the MARTINI 
model considers only reversible equilibrium behaviors. For 
example, another feature of a fully hydrated phosphatidyl¬ 
choline bilayer, observed in cooling protocols of multilamellar 
vesicles and planar bilayers, is the appearance of the so-called 
ripple phase, where the surface of the bilayer appears corru¬ 
gated [2]. The ripple structure is not observed in the equilib¬ 
rium behavior of the MARTINI model, but it is worth study¬ 
ing whether this structure and others emerge out of equilib¬ 
rium, as metastable states. Moreover, both tilt as well as the 
rippled structure disappear upon addition of a small percent¬ 
age of impurities such as cholesterol [42,43], which is another 
behavior worthy of future study. 

Concerning our finding of first-order character in the liq- 
uid-hexatic transition, experiments have found such behavior 
in five-layer thick smectic liquid-crystal films [44]. Further, re¬ 
cent simulations indicate that the hexatic to liquid transition 
is in fact a first-order transition in two-dimensional systems 
of hard disks [35]. Other work on two-dimensional systems 
shows that the order of the hexatic to liquid transition de¬ 
pends on the steepness of the repulsive interactions between 
particles [45]. Specifically, it was found that for soft disks 
with repulsive potentials the transition is continuous for 
n < 6 while it is first order for n > 6 . Therefore, the order of 
the hexatic-to-liquid phase transition is a system-dependent 
property. Indeed, a different model for lipid bilayers with soft 
repulsive interactions [41,46] was found by Rodgers et al. to 
change continuously between ordered and disordered struc¬ 
tures. In that model, it is not even clear that a transition 
exists between distinct phases. 

This system dependence is more than a theoretical cu¬ 
riosity. Experiments can vary the presence and order of a 
transition by varying membrane composition. For composi¬ 
tions where the transition remains first order, pre-melting of 
the ordered phase and related phenomena become relevant. 
In particular, large enough membrane-bound solutes, (e.g., 
transmembrane proteins) can form surfaces that disfavor the 
ordered phase even at thermodynamic conditions where the 
ordered phase is stable. At such conditions, a microscopic pre¬ 
melting layer with disordered membrane structure will then 
necessarily surround these solutes. Several such neighboring 
solutes will then assemble to minimize line tension. 

This mechanism of self assembly is the subject of the sub¬ 
sequent paper [1]. The strength of the assembling force will 




Fig. 4. Hexatic phase; (A) Dislocations (5-7 disclination pairs) in a 25X25 nm^ patch of bilayer taken from a lOTXllOnm^ sample in the ordered phase at 294 K. 
The tail-end particles of one monolayer are shown in gray. Particles with seven neighbors are highlighted in red, while particles with five neighbors are highlighted in blue. 
There exist several unbound dislocations as expected for a hexatic phase. These dislocations appear to diffuse freely throughout the system as seen from three different 
configurations separated by 600 ps. (B) One-dimensional pair correlation function g{x) showing decay faster than . (C) Orientational correlation function Gq{x) 

showing quasi-long-range correlations at T = 294 K and short-range correlations at T = 312 K. Values of Gq{x) smaller than 10“^ amount to statistical noise for the 
system sizes considered. (D) Enlargement of Panel C to show the slow decay of the quasi-long-range orientational correlations at T = 294 K. 


depend upon the temperature, lateral pressure and the com¬ 
position of the membrane, all of which affect the nature of 
the order-disorder transition. Without its first-order charac¬ 
ter, solute-induced pre-melting layers will not exist, line ten¬ 
sion will not exist, and this force of assembly will vanish. We 
think the possibility of moving between regimes with differ¬ 
ing transition order implies an important degree of complexity 
and correlated behavior yet to be explored. 

Materials and Methods 

We use the MARTINI coarse-grained force field [4] to model a lipid bilayer system in 
which four carbon atoms (or equivalent) are approximated as one coarse-grained bead. 
This model has an explicit solvent with approximately four water molecules scaled to 
one solvent particle. Dipalmitoyl phosphatidylcholine (DPPC) is chosen as the model 
lipid species. The MARTINI model uses the Lennard-Jones potential for non-bonded 
interactions. The cut-off for these interactions is 1.2 nm. The GROMACS shifting 
function [47] is used in the range 0.99-1.2 nm. Bond and angle energies are mod- 

— 1 —2 

eled as harmonic potentials with associated force constants 1250 kJmol nm 
— 1 —2 

and 25kJmol rad respectively. We use a time step of 30fs, within the rec¬ 
ommended range of 20-40fs for this model. Simulations are performed using the 
GROMACS molecular dynamics package [48] in the ensemble with fixed numbers of 
particles, temperature, and stress tensor components for the system [4,49,50]. Ther¬ 
mostats and barostats were used to control temperature and pressure, and checks 
were performed to assure that different thermostats and barostats yielded similar re¬ 
sults [51]. The compressibility is set to 3X10 ^ bar A lateral pressure of zero 
is maintained by choosing the diagonal components of the stress tensor to be 1 bar. 
The electrostatic interactions are shifted to zero in the range 0-1.2 nm. A dielectric 
constant of 15 is used for screening of electrostatic interactions. 


The flat interface is stabilized by juxtaposing an ordered bilayer equilibrated at 
285 K and zero lateral pressure with a disordered bilayer equilibrated at the same 
conditions corresponding to the cooling and heating curves of the hysteresis loop in 
Fig. 2, respectively. The system thus constructed is equilibrated in the ensemble 
with fixed temperature, volume and numbers of particles. This ensemble allows for 
maintaining an area per lipid intermediate between the two phases, thus stabilizing 
the interface. 

The translational and orientational correlation functions for the ordered phase are 
calculated for a system containing 51,200 lipids solvated by 800,000 water molecules 
(1.4 million particles). This system is prepared by first equilibrating a smaller sys¬ 
tem containing 3,200 lipids to make a hydrated bilayer membrane of area roughly 
2 

26X27nm .A configuration from this equilibrated system is then periodically repli¬ 
cated to produce an initial configuration for the simulated system, which has a mem- 
2 

brane area of roughly 107X110nm . The system is then relaxed for about 11 mi¬ 
croseconds resulting in a hexatic phase. 

To identify 5-7 disclination pairs, we use Shewchuk’s code. Triangle [52,53]. 
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Supporting Information 

for is a first-order hexatic to liquid phase transition” 

by Shachi Katira, Kranthi K. Mandadapu, Suriyanarayanan Vaikuntanathan, Berend 
Smit and David Chandler 


In this Supporting Information, we provide evidence that the different planes containing the different beads along the length 
of the chain also have short-range translational order and quasi-long-range orientational order, thus indicating that the entire 
hydrophobic part of the bilayer is in the hexatic phase. We also discuss the relaxation of the bilayer system into the hexatic 
phase and the corresponding relaxation of the translational, pair and orientational correlation functions for the million-particle 
system. 

Along the lipid chain 

In addition to the tail-end particles (C4), we calculate pair correlation functions for non-tail-end particles C2 and C3 on the 
lipid chains. The planes containing these particles show fewer unbound dislocations compared to the tail-end particles C4 as 
shown in Fig. SI. The amount of disorder therefore decreases higher along the lipid chains, consistent with the disorder gradient 
observed in experiments [37] and explained by theory [38]. Despite the gradient in disorder, the pair correlation functions, 
g{x)^ seen in Fig. S2A appear to decay exponentially for all three planes of particles. The orientational correlation functions 
also exhibit quasi-long-range behavior Fig. S2B. Moreover, the dislocations in the C2 and C3 planes also diffuse freely in the 
time scale of the simulation, albeit slower than those in the C4 plane. These data show that all three planes of particles are 
hexatic. 
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Fig. SI. 5-7 disclination pairs for planes containing C2 (left), C3 (center) and C4 (right) particles in a 107X110 nm^ bilayer in the ordered phase at T=294 K. The tail-end 
particles of one monolayer are shown in gray. Enlarged red and blue circles render particles with seven and five neighbors, respectively. For each plane, there exist unbound 
dislocations as expected for a hexatic phase, and the number of dislocations increases towards the tail-end. The box drawn in for the C4 panel shows the region illustrated in 
Panel A of Fig. 4 in the main text. 


Relaxation 

In Fig. S3, we show the time evolution of the pair correlation function for the most ordered C2 and the least ordered C4 
planes. The correlations observed at 0.5 microseconds reflect the preparation of the system by periodically replicating a small 
ordered patch several times in the x and y directions. Specifically, notice the broad peaks at those early times, as the initial 
configuration is a lattice, effectively a crystal, with large unit cells. As time progresses, the lattice dissipates, positions become 
increasingly uncorrelated, and exponential decay emerges for g{x) — 1 in both C2 and C4 planes. 

In Fig. S4, we show the time evolution of the orientational correlation function G6 (t) for the C2 and C4 planes of particles. 
The long-range correlations introduced by periodically replicating the system gradually diminish and the system shows a 
quasi-long-range behavior at 6.3 /xs. 
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Fig. S3. Translational pair correlation function during relaxation to the hexatic phase, starting from initial non-equilibrium state. (A) g(x) — 1 for tail particles C4 at three 
different times. (B) g{x) — 1 for tail particles C2 at three different times. 
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Fig. S4. Orientational pair correlation function during relaxation to the hexatic phase, starting from initial non-equilibrium state. (A) Gq{x) for tail particles C4 at three 
different times. (B) Gq{x) for tail particles C2 at three different times. 
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